(library (grid-gen symath)

  (export solve-eq
          solve-sys
          diff-eq
          jakobi-matrix
          simplify-form
          sort-form
          get-all-vars
          comp-numbers
          expand-form
          ;;infix-out
          ;;get-right-ev
          ;;get-left-ev
          easy-simplify-form
          simp-ex
          simp-ex2
          simplify4-form
          norm-form
          sym-det
          sym-vec-mul
          sym-vec-mod
    ;;      grad ;; untested
          expand-to-wires
          subst-vars
          l->s
  )

;;(use-modules (statprof))
(import (guile)(rnrs);;(chezscheme)
        (grid-gen main)
        (grid-gen vars)
        (grid-gen code)
        (grid-gen math)
        ;;(except (grid-gen symath) (grid-gen math))
        (except (srfi srfi-1) map)
        ;;(srfi :28)
)

;;  $X -- any form
;;  %X -- form contains variables only
;; %!X -- variable
;; %%X -- form with a variable
;; !%X -- form without any variable
;;  &X -- number
;; +&X -- positive number
;; -&X -- negative &X
;; &&X -- form contains numbers only
;; !&X -- not a number
;; _N  -- number N

(define expand-templates
  '(
    ((inv (inv $1)) $1) ;;+
    ((neg (neg $1)) $1) ;;+

    ((- $1 $2) (+ $1 (neg $2))) ;;+
    ((/ $1 $2) (* $1 (inv $2))) ;;+
    
    ((* $1 _-1) (neg $1)) ;;+
    ((* _-1 $1) (neg $1)) ;;+

    ((neg (* !&1 $2)) (* (neg !&1) $2)) ;;+
    ((neg (* $2 !&1)) (* (neg !&1) $2)) ;;+
    ((neg (+ $1 $2)) (+ (neg $1) (neg $2))) ;;+

    ((inv (* $1 $2)) (* (inv $1) (inv $2))) ;;+
    
    ((* $1 (+ $2 $3)) (+ (* $1 $2) (* $1 $3))) ;;+
    ((* (+ $1 $2) $3) (+ (* $3 $1) (* $3 $2))) ;;+
   )
)

(define simplify1-templates
  '(
    ((* $1 _1) $1)  ;;+
    ((* $1 _0) 0.0) ;;+
    ((+ $1 _0) $1)  ;;+
    ((* _1 $1) $1)  ;;+
    ((* _0 $1) 0.0) ;;+
    ((+ _0 $1) $1)  ;;+
    ((** $1 _1) $1)    ;;+
    ((** $1 _0) 1.0)   ;;+
    ((** $1 _-1) (inv $1))  ;;+

    ((sqrt $1) (** $1 0.5)) ;;+
    ((inv (inv $1)) $1)     ;;+
    ((neg (neg $1)) $1)
    ((- $1 $2) (+ $1 (neg $2))) ;;+
    ((/ $1 $2) (* $1 (inv $2))) ;;+
    ;; ((neg $1) (* -1.0 $1))      ;;+
    
    ((* &1 !&2) (* !&2 &1)) ;;+
    ((+ &1 !&2) (+ !&2 &1)) ;;+
    ((* (* !&1 &2) (* !&3 &4)) (* (* !&1 !&3) (* &2 &4))) ;;+

    ((* (* $1 &2) $3) (* $1 (* $3 &2))) ;;+
    ((+ (+ $1 &2) &3) (+ $1 (+ &2 &3))) ;;+

    ((+ (* $1 +&2) (* $3 -&2)) (* +&2 (+ $1 (neg $3)))) ;;+
    ((+ (* $1 -&2) (* $3 +&2)) (* +&2 (+ $3 (neg $1)))) ;;+
    
    ((** (neg $1) _2) (** $1 2.0)) ;;+
    
    ((* (neg $1) $1) (neg (** $1 2.0)))
    ((* $1 (neg $1)) (neg (** $1 2.0)))
   )
)

(define simplify2-templates
  '(
    ((+ $1 $1) (* 2.0 $1))  ;;+
    ((* $1 (inv $1)) 1.0)   ;;+
    ((* $1 $1) (** $1 2.0)) ;;+
    ((* $1 (* $2 (inv $1))) $2) ;;+
    ((* $1 (* (inv $1) $2)) $2) ;;+
    ((* (inv $1) (* $2 $1)) $2) ;;+
    ((* (inv $1) (* $1 $2)) $2) ;;+
    ((* (* $2 $1) (inv $1)) $2) ;;+
    ((* (* $1 $2) (inv $1)) $2) ;;+
    ((* $1 (inv (* $1 $2))) (inv $2)) ;;+
    
    ((inv (* $1 (inv $2))) (* $2 (inv $1)))         ;;+
    ((+ (* $1 $2) $1) (* $1 (+ $2 1.0)))            ;;+
    ((+ (* $1 $2) (* $3 $2)) (* $2 (+ $1 $3)))      ;;+
    ((* (* $1 $2) (inv (* $1 $3))) (* $2 (inv $3))) ;;+
    ((* (** $1 $2) (inv (* $1 $3))) (* (** $1 (+ $2 -1.0)) (inv $3))) ;;+
    ((+ (* $1 $2) (* $3 (* $4 $2))) (* $2 (+ $1 (* $3 $4))))          ;;+
    ((* (* (** $1 $2) $4) (inv (* $1 $3))) 
     (* (* (** $1 (+ $2 -1.0)) $4) (inv $3)))         ;;+
    ((* (* (** $1 $2) $4) (inv (* (** $1 $5) $3))) 
     (* (* (** $1 (+ $2 (neg $5))) $4) (inv $3)))     ;;+
    ((+ (* $3 (** $1 $2)) (* $5 (* $1 $4)))
     (* $1 (+ (* (** $1 (+ $2 -1.0)) $3) (* $4 $5)))) ;;+
    ((+ $1 (* $1 $2)) (* $1 (+ $2 1.0)))       ;;+

    ((* (** $1 $2) $1) (** $1 (+ $2 1.0)))     ;;+
    ((** (** $1 $2) $3) (** $1 (* $2 $3)))     ;;+
    ((+ (* $1 $2) (* $1 $3)) (* $1 (+ $2 $3))) ;;+
    ((+ (* $1 $2) (* $3 $2)) (* (+ $1 $3) $2)) ;;+
    ((* (** $1 $2) (** $1 $3)) (** $1 (+ $2 $3)))        ;;+
    ((* $1 (* $2 (* $1 $3))) (* (** $1 2.0) (* $2 $3)))  ;;+
    ((* $1 (* (** $1 $2) $3)) (* (** $1 (+ $2 1.0)) $3)) ;;+
    ((* (** $1 $2) (* $1 $3)) (* (** $1 (+ $2 1.0)) $3)) ;;+
    ((* (** $1 $2) (* (** $1 $3) $4)) (* (** $1 (+ $2 $3)) $4))        ;;+
    ((* $1 (* $2 (* (** $1 $3) $4))) (* $2 (* (** $1 (+ $3 1.0)) $4))) ;;+

    ((+ (+ (* (+ $1 $2) $3) $1) $4)
     (+ (* $1 (+ $3 1.0)) (+ (* $2 $3) $4)))
    ((+ (* $1 $2) (* (** $1 $3) $4))
     (* $1 (+ $2 (* (** $1 (+ $3 -1.0)) $4))))
    ((* (** $1 $2) (* (** (* $1 $3) $4) $5))
     (* (** $1 (+ $2 $4)) (* (** $3 $4) $5)))
    ((* (** $1 $2) (* $3 (* $1 $4))) (* (** $1 (+ $2 1.0)) (* $3 $4)))
    ((** (* $1 (** $2 $3)) $4) (* (** $1 $4) (** $2 (* $3 $4))))
    ((inv (* (** $1 $2) $3)) (* (** $1 (neg $2)) (inv $3)))
    ((inv (** $1 $2)) (** $1 (neg $2)))                    
    ((+ $1 (** $1 $2)) (* $1 (+ (** $1 (+ $2 -1.0)) 1.0)))
    ((* (+ $1 (* $2 (* (** $3 $4) $5))) (** $3 $6))
     (+ (* $1 (** $3 $6)) (* $2 (* (** $3 (+ $4 $6)) $5))))    
    ((** (* $1 $2) $3) (* (** $1 $3) (** $2 $3)))

    ((+ (* (** $1 $2) $3) (* $1 $4))
     (* $1 (+ (* (** $1 (+ $2 -1.0)) $3) $4)))
    ((+ (* (** $1 $2) $3) (+ (* $1 $4) $5))
     (+ (* $1 (+ (* (** $1 (+ $2 -1.0)) $3) $4)) $5))
    ((+ (* (+ $1 $2) $3) (+ $1 $4)) (+ (* $1 (+ $3 1.0)) (+ (* $2 $3) $4)))
    ((+ (* (* $1 $2) $3) (* $2 $4)) (* $2 (+ (* $1 $3) $4)))
    ((* (* (* $1 $2) $3) (inv (* $4 $1))) (* (* $2 $3) (inv $4)))
    ((+ (* $1 $2) (+ $3 (* $1 $4))) (+ (* $1 (+ $2 $4)) $3))    
    ((+ (* $1 (* $2 $3)) (* $2 $4)) (* $2 (+ (* $1 $3) $4)))
    
;;;;;;;;;;;;;;;;;
    ((* (inv $1) $2) (* $2 (inv $1)))
    ((** (inv $1) $2) (inv (** $1 $2)))
    ((* $1 (* $2 (inv $3))) (* (* $1 $2) (inv $3)))
    ((* (* $1 (inv $2)) (inv $3)) (* $1 (inv (* $2 $3))))
    
;;;;;;;;;;;;;;;;;
    ((+ (* $1 (inv $2)) (* $3 (inv $4))) 
     (* (+ (* $1 $4) (* $3 $2)) (inv (* $2 $4))))
    ((** (* $1 (inv $2)) $3) (* (** $1 $3) (inv (** $2 $3))))
    ((+ $1 (inv $2)) (* (+ (* $1 $2) 1.0) (inv $2)))
    ((+ $1 (* $2 (inv $3))) (* (+ (* $1 $3) $2) (inv $3)))
    ((+ (* $1 (inv $2)) $3) (* (+ $1 (* $3 $2)) (inv $2)))
    
    ((+ (* $1 (* $2 $3)) (* $4 (* $5 $3))) (* (+ (* $1 $2) (* $4 $5)) $3))
    ((+ (* (+ $1 (neg $2)) $3) (* (+ (neg $1) $2) $4))
     (* (+ $1 (neg $2)) (+ $3 (neg $4))))
    
    ((* (inv !&1) (inv !&2)) (inv (* !&1 !&2)))

    ((neg (* $1 &2)) (* $1 (neg &2)))
    ((neg (* &2 $1)) (* $1 (neg &2)))
        
    ((+ (+ !&1 &2) !&3) (+ (+ !&1 !&3) &2))
    ((* !&1 (* !&2 &3)) (* (* !&1 !&2) &3))

    ((inv (* $1 &2)) (* (inv $1) (inv &2)))
    
    ((+ (* $1 $2) (neg (* $3 $2))) (* $2 (+ $1 (neg $3))))
    ((+ (* (+ (neg $1) $2) $3) (* (+ $1 (neg $2)) $4))           
     (* (+ (neg $1) $2) (+ $3 (neg $4))))

    ((* (* (+ (neg $1) $2) $3) (inv (* (+ $1 (neg $2)) $4)))
     (neg (* $3 (inv $4))))
    ((* (* (+ $2 (neg $1)) $3) (inv (* (+ $1 (neg $2)) $4)))
     (neg (* $3 (inv $4))))
    ((* (* (+ (neg $1) $2) $3) (inv (* (+ (neg $2) $1) $4)))
     (neg (* $3 (inv $4))))
    ((* (* (+ $2 (neg $1)) $3) (inv (* (+ (neg $2) $1) $4)))
     (neg (* $3 (inv $4))))

    ((* (* (+ $1 (neg $2)) $3) (inv (* (+ (neg $1) $2) $4)))
     (neg (* $3 (inv $4))))
    ((* (* (+ (neg $2) $1) $3) (inv (* (+ (neg $1) $2) $4)))
     (neg (* $3 (inv $4))))
    ((* (* (+ $1 (neg $2)) $3) (inv (* (+ $2 (neg $1)) $4)))
     (neg (* $3 (inv $4))))
    ((* (* (+ (neg $2) $1) $3) (inv (* (+ $2 (neg $1)) $4)))
     (neg (* $3 (inv $4))))

    ((* (+ (neg $1) $2) (inv (+ $1 (neg $2)))) -1.0)
    ((* (+ $2 (neg $1)) (inv (+ $1 (neg $2)))) -1.0)
    ((* (+ $1 (neg $2)) (inv (+ (neg $1) $2))) -1.0)
    ((* (+ (neg $2) $1) (inv (+ (neg $1) $2))) -1.0)

    ((+ (* $1 (+ (* $2 $3) $4)) (* $2 $5)g)                     
     (+ (* $2 (+ (* $1 $3) $5)) (* $1 $4)))
     
    ((* (* (* $2 $1) $3) (* $4 (* $5 $1)) $6)
     (* (** $1 2.0) (* (* $2 $3) (* (* $4 $5) $6))))
     
    ((+ (* $1 $2) (neg (* $2 $3))) (* $2 (+ $1 (neg $3))))
    
    ((* (* $1 $2) (inv (* $3 $1))) (* $2 (inv $3)))
    ((+ (* $1 $2) (neg (* $1 $3))) (* $1 (+ $2 (neg $3))))
    ((neg (+ (neg $1) $2)) (+ $1 (neg $2)))
    ((+ (* (+ (* $1 $2) $3) $4) (neg (** $1 $5)))
     (+ (* $1 (+ (* $2 $4) (** $1 (+ $5 -1.0)))) (* $3 $4)))
   )
)

(define sort-templates
  '(
    ((+ !%1 %%1) (+ %%1 !%1))
    ((* !%1 %%1) (* %%1 !%1))

    ((* (neg %%1) !%2) (* %%1 (neg !%2)))

    ((+ !%1 (+ %%1 $2)) (+ %%1 (+ !%1 $2)))
    ((+ (+ %%1 $1) $2)  (+ %%1 (+ $1 $2)))
    
    ((+ (* %%1 $2) (+ (* %%1 $3) $4)) (+ (* %%1 (+ $2 $3)) $4))
    
    ((+ (* %%1 $1) (* %%1 $2)) (* %%1 (+ $1 $2)))

    ((* !%1 (* %%1 $2)) (* %%1 (* !%1 $2)))
    ((* (* %%1 $1) $2) (* %%1 (* $1 $2)))
    
    ((* %%1 %%1) (** %%1 2.0))
    ((* %%1 (* %%1 $2)) (* (** %%1 2.0) $2))
    ((* %%1 (** %%1 $2)) (** %%1 (+ $2 1.0)))
    
    ((* &1 !&1) (* !&1 &1))
    ((* &1 (* !&2 !&3)) (* !&2 (* !&3 &1)))
    
    ((+ &1 !&1) (+ !&1 &1))
    ((+ !&1 (+ !&2 &3)) (+ (+ !&1 !&2) &3))

   )
)

(define simp1-templates
  '(
    ((* $1 _1) $1)  ;;+
    ((* $1 _0) 0.0) ;;+
    ((+ $1 _0) $1)  ;;+

    ((* _1 $1) $1)  ;;+
    ((* _0 $1) 0.0) ;;+
    ((+ _0 $1) $1)  ;;+

    ((** $1 _1) $1)    ;;+
    ((** $1 _0) 1.0)   ;;+
    ((** $1 _-1) (inv $1))  ;;+
    
    ((- $1 $1) 0.0)
    ((+ $1 (neg $1)) 0.0)
    ((+ (neg $1) $1) 0.0)

    ((/ $1 $1) 1.0)
    ((* $1 (inv $1)) 1.0)
    ((* (inv $1) $1) 1.0)

    ((sqrt $1) (** $1 0.5)) ;;+
    ((inv (inv $1)) $1)     ;;+
    ((- $1 $2) (+ $1 (neg $2))) ;;+
    ((/ $1 $2) (* $1 (inv $2))) ;;+
    ((neg $1) (* -1.0 $1))      ;;+
   )
)


(define compact-templates
  '(
    ((+ (* $1 $2) (* $1 $3)) (* $1 (+ $2 $3)))
    ((* $1 (neg !&2)) (neg (* $1 !&2)))
    ((* (neg !&1) $2) (neg (* !&1 $2)))
    ((inv (neg $1)) (neg (inv $1)))
    ((* (* $1 (inv $2)) (* $3 (inv $4))) (* (* $1 $3) (inv (* $2 $4))))
    ((+ (* $1 (inv $2)) (* $3 (inv $4))) 
     (* (+ (* $1 $4) (* $3 $2)) (inv (* $2 $4)))
    )
   )
)

(define solve-templates
  '(
    ((+ %!1 !%1) (neg !%1))
    ((+ !%1 %!1) (neg !%1))

;;    ((+ (* !%1 %!1) !%2) (* (neg !%2) (inv !%1)))
    ((+ (* %!1 !%1) !%2) (* (neg !%2) (inv !%1)))
    
    ((* %!1 !%1) 0.0)
    
    ((+ (** %!1 !%2) !%1) (** (neg !%1) (neg !%2)))
    ((+ (* (** %!1 !%1) !%2) !%3) (** (* (neg !%3) (inv !%2)) (neg !%1)))
    
    ((+ (* !%1 (inv (** %!1 !%2))) !%3)
     (** (* !%1 (inv (neg !%3))) (inv !%2))
    )

    ((+ (* (inv (** %!1 !%2)) !%1) !%3)
     (** (* !%1 (inv (neg !%3))) (inv !%2))
    )

    ((+ !%3 (* !%1 (inv (** %!1 !%2))))
     (** (* !%1 (inv (neg !%3))) (inv !%2))
    )

    ((+ !%3 (* (inv (** %!1 !%2)) !%1))
     (** (* !%1 (inv (neg !%3))) (inv !%2))
    )

;;    ((+ !%1 (** %!1 !%2)) (** (neg !%1) (neg !%2)))
;;    ((+ !%2 (* %!1 !%1)) (* (neg !%2) (inv !%1)))
;;    ((+ !%2 (* !%1 %!1)) (* (neg !%2) (inv !%1)))
;;    ((+ !%3 (* (** %!1 !%1) !%2)) (** (* (neg !%3) (inv !%2)) (neg !%1)))
;;    ((+ (* !%2 (** %!1 !%1)) !%3) (** (* (neg !%3) (inv !%2)) (neg !%1)))
;;    ((+ !%3 (* !%2 (** %!1 !%1))) (** (* (neg !%3) (inv !%2)) (neg !%1)))
   )
)

(define math-num-forms
  '(
    (+ &1 &2)
    (- &1 &2)
    (* &1 &2)
    ;;(/ &1 &2)
    (** &1 &2)
    (sin &1)
    (cos &1)
    (tan &1)
    (asin &1)
    (acos &1)
    (atan &1)
    (log &1)
    (exp &1)
    (sqrt &1)
    (abs &1)
    (inv &1)
    (neg &1)
    (abs &1)
    (sign &1)
   )
)

(define norm1-templates
  '(
    ((** $1 _-1) (inv $1))

    ((+ $1 _-1) (- $1 1.0))
    ((+ _-1 $1) (- $1 1.0))

    ((* $1 (inv $2)) (/ $1 $2))
    ((* (inv $2) $1) (/ $1 $2))
    ((* $1 (/ _1 $2)) (/ $1 $2))
    ((* (/ _1 $2) $1) (/ $1 $2))

    ((+ $1 (neg $2)) (- $1 $2))
    ((+ (neg $2) $1) (- $1 $2))

    ((* $1 _-1) (* -1.0 $1))
    ((* (* _-1 $1) $2) (* -1.0 (* $1 $2)))
    ((* $1 (* _-1 $2)) (* -1.0 (* $1 $2)))

    ((+ $1 (* _-1 $2)) (- $1 $2))
    ((+ (* _-1 $2) $1) (- $1 $2))

    ((* _-1 (- $1 $2)) (- $2 $1))
    
    ((** $1 _0.5) (sqrt $1))

    ((neg $1) (* -1.0 $1))
    ((inv $1) (/ 1.0 $1))
    
    ((/ $1 _2) (* 0.5 $1))
    
   )
)

(define norm2-templates
  '(
    ((** $1 _2) (* $1 $1))
    ((+ (+ $1 $2) $3) (+ $1 $2 $3))
    ((+ $1 (+ $2 $3)) (+ $1 $2 $3))
    ((* (* $1 $2) $3) (* $1 $2 $3))
    ((* $1 (* $2 $3)) (* $1 $2 $3))
   )
)

(define diff-templates
  '(
    ((diff %!1) 1.0)
    ((diff !%1) 0.0)

    ((diff (inv $1)) (neg (* (diff $1) (inv (** $1 2.0)))))
    ((diff (neg $1)) (neg (diff $1)))

;;    ((diff (** $1 _2)) (* 2.0 $1))
    ((diff (** $1 &&2)) (* (* &&2 (** $1 (+ &&2 -1.0))) (diff $1)))
    
;;    ((diff (** &&1 %%1)) (* (* (** &&1 %%1) (log &&1)) (diff %%1)))

    ((diff (+ $1 $2)) (+ (diff $1) (diff $2)))
    ((diff (* $1 $2)) (+ (* (diff $1) $2) (* $1 (diff $2))))

    ((diff (exp $1)) (* (exp $1) (diff $1)))
    ((diff (log $1)) (* (inv $1) (diff $1)))
    ((diff (sin $1)) (* (cos $1) (diff $1)))
    ((diff (cos $1)) (* (neg (sin $1)) (diff $1)))
    ((diff (tan $1)) (* (inv (** (cos $1) 2.0)) (diff $1)))
    ((diff (asin $1)) (* (inv (sqrt (+ 1.0 (neg (** $1 2.0))))) (diff $1)))
    ((diff (acos $1)) 
     (* (neg (inv (sqrt (+ 1.0 (neg (** $1 2.0)))))) (diff $1))
    )
    ((diff (atan $1)) (* (inv (+ 1.0 (** $1 2.0))) (diff $1))
    )
    
    ((diff (sqrt $1)) (diff (** $1 0.5)))
    
    ((diff (** $1 !&2))
     (* (** $1 !&2)
        (+ (* (diff !&2) (log $1))
           (* (* !&2 (diff $1))
              (inv $1)
           )
        )
     )
    )
  )
)


(define chop-10 (lambda (l)
  l
;;  (list-head l (min 5 (length l)))
))


(define is-tmpl? (lambda (s)
  (and (symbol? s)
       (string-prefix? "$" (sym->str s))
  )
))

(define is-proc-tmpl? (lambda (s)
  (and (symbol? s)
       (string-prefix? "%" (sym->str s))
  )
))

(define is-proc2-tmpl? (lambda (s)
  (and (symbol? s)
       (string-prefix? "%%" (sym->str s))
  )
))

(define is-proc3-tmpl? (lambda (s)
  (and (symbol? s)
       (string-prefix? "!%" (sym->str s))
  )
))

(define is-proc4-tmpl? (lambda (s)
  (and (symbol? s)
       (string-prefix? "&" (sym->str s))
  )
))


(define is-proc5-tmpl? (lambda (s)
  (and (symbol? s)
       (string-prefix? "_" (sym->str s))
  )
))

(define is-proc6-tmpl? (lambda (s)
  (and (symbol? s)
       (string-prefix? "!&" (sym->str s))
  )
))

(define is-proc7-tmpl? (lambda (s)
  (and (symbol? s)
       (string-prefix? "%!" (sym->str s))
  )
))

(define is-proc8-tmpl? (lambda (s)
  (and (symbol? s)
       (string-prefix? "&&" (sym->str s))
  )
))

(define is-proc-minus? (lambda (s)
  (and (symbol? s)
       (string-prefix? "-&" (sym->str s))
  )
))

(define is-proc-plus? (lambda (s)
  (and (symbol? s)
       (string-prefix? "-&" (sym->str s))
  )
))

(define is-proc-plus? (lambda (s)
  (and (symbol? s)
       (string-prefix? "+&" (sym->str s))
  )
))


(define get-proc5-num (lambda (s)
  (let* ((ss (sym->str s))
         (ssl (string-length ss))
        )
    (string->number (substring ss 1 ssl))
  )
))

(define cont-only (lambda (v l)
  (or (and (or (symbol? v) (procedure? v)) (in-list? v l))
      (and (list? v)
           (eq? (length v) 2)
           (or (eq? (car v) 'neg) (eq? (car v) 'inv))
           (cont-only (cadr v) l)
      )
  )
))

(define cont-only-nums (lambda (v)
  (letrec ((c-o-n-l (lambda (l)
                      (or (null? l)
                          (and (cont-only-nums (car l))
                               (c-o-n-l (cdr l))
                          )
                      )
                    )
           )
          )
    (or (number? v)
        (and (list? v)
             (not (null? v))
             (c-o-n-l (cdr v))
        )
    )
  )
))

(define cont-one (lambda (v l)
  (cond ((or (symbol? v) (procedure? v))
         (in-list? v l)
        )
        ((and (list? v) (not (null? v)))
         (or (cont-one (car v) l)
             (cont-one (cdr v) l)
         )
        )
        (#t #f)
  )
))

(define form-vars (lambda (f t vl vrs)
  (let ((rrr
  (let ((tv (hash-ref vrs t #f))
        (ret (lambda ()
          (hash-set! vrs t f)
          vrs
        ))
       )
    (cond (tv
           (if (equal? f tv)
               vrs
           )
          )
          ((and (is-proc-minus? t) (number? f) (< f 0.0))
           (let* ((ss (sym->str t))
                  (ss1 (string-take-right ss (- (string-length ss) 1)))
                  (ss2 (string-append "+" ss1))
                  (ts2 (string->symbol ss2))
                  (tv2 (hash-ref vrs ts2 #f))
                 )
             (if tv2
                 (if (= f (* -1.0 tv2))
                     (ret)
                 )
                 (ret)
             )
           )
          )
          ((and (is-proc-plus? t) (number? f) (>= f 0.0))
           (let* ((ss (sym->str t))
                  (ss1 (string-take-right ss (- (string-length ss) 1)))
                  (ss2 (string-append "-" ss1))
                  (ts2 (string->symbol ss2))
                  (tv2 (hash-ref vrs ts2 #f))
                 )
             (if tv2
                 (if (= f (* -1.0 tv2))
                     (ret)
                 )
                 (ret)
             )
           )
          )
          ((is-proc8-tmpl? t)
           (if (cont-only-nums f)
               (ret)
           )
          )
          ((is-proc7-tmpl? t)
           (if (in-list? f vl)
               (ret)
           )
          )
          ((is-proc6-tmpl? t)
           (if (not (number? f))
               (ret)
           )
          )
          ((is-proc5-tmpl? t)
           (if (and (number? f)
                    (= (string->number (number->string f))
                       (get-proc5-num t)
                    )
               )
               vrs
           )
          )
          ((is-proc4-tmpl? t)
           (if (number? f)
               (ret)
           )
          )
          ((is-proc3-tmpl? t)
           (if (not (cont-one f vl))
               (ret)
           )
          )
          ((is-proc2-tmpl? t)
           (if (cont-one f vl)
               (ret)
           )
          )
          ((is-proc-tmpl? t)
           (if (cont-only f vl)
               (ret)
           )
          )
          ((is-tmpl? t) (ret))
          ((equal? f t) 
           vrs
          )
          ((and (list? f) (not (null? f))
                (list? t) (not (null? t))
           )
           (let* ((fv1 (form-vars (car f) (car t) vl vrs))
                  (fv2 (form-vars (cdr f) (cdr t) vl vrs))
                 )
         ;;     (format #t "F: ~a\nT: ~a\nFV1: ~a\nFV2: ~a\n" f t fv1 fv2)
              (if (and (not (unspecified? fv1))
                       (not (unspecified? fv2))
                  )
                  vrs
              )
           )
          )
    )
  )))
;;    (format #t "XF: ~a\nXT: ~a\nXR: ~a\n" f t rrr)
    rrr
  )
))


(define subst-vars (lambda (t v)
  (cond ((or (symbol? t) (procedure? t)) ;;;;;;;;;;;;;;;;;;;;;;;;;
         (hash-ref v t t)
        )
        ((and (list? t) (not (null? t)))
         (cons (subst-vars (car t) v)
               (subst-vars (cdr t) v)
         )
        )
        (#t t)
  )
))

(define replace-simple-form (lambda (f t vl)
  (let* ((tv (form-vars f (car t) vl (make-hash-table)))
        )
  ;;  (format #t "RSF: ~a\nRST: ~a\nRSV: ~a\n" f (car t) tv)
    (if (unspecified? tv)
        (cond ((symbol? f) f)
              ((procedure? f) f)
              ((number? f) f)
              ((null? f) f)
              ((unspecified? f) f)
              ((list? f)
               (cons (replace-simple-form (car f) t vl)
                     (replace-simple-form (cdr f) t vl)
               )
              )
              (#t (int-err "Cannot expand form value"))
        )
        (let ((rrr
        (subst-vars (cadr t) tv)
        ))
    ;;    (format #t "RES: ~a\n\n" rrr)
        rrr
        )
    )
  )
))

(define replace-simple-forms (lambda (f tl vl)
  (if (null? tl)
      f
      (replace-simple-forms
        (replace-simple-form f (car tl) vl)
        (cdr tl)
        vl
      )
  )
))

(define expand-to-wire+ (lambda (v)
  (if (eq? (length (cdr v)) 2)
      v
      (list (car v)
            (cadr v)
            (expand-to-wire+ (cons (car v) (cddr v)))
      )
  )
))

(define expand-to-wire- (lambda (v)
  (if (eq? (length (cdr v)) 2)
      v
      (list (car v)
            (expand-to-wire- 
              (cons (car v)
                    (reverse (cdr (reverse (cdr v))))
              )
            )
            (car (reverse v))
      )
  )
))

(define expand-to-wires (lambda (v)
  (cond ((and (list? v) (not (null? v))) 
         (cond ((or (eq? (car v) '+)
                    (eq? (car v) '*)
                )
                (expand-to-wire+ (cons (car v) (map expand-to-wires (cdr v))))
               )
               ((or (eq? (car v) '-)
                    (eq? (car v) '/)
                )
                (expand-to-wire- (cons (car v) (map expand-to-wires (cdr v)))) 
               )
               (#t (cons (car v) (map expand-to-wires (cdr v))))
         )
        )
        (#t v)
  )
))


(define replace-form (lambda (f t vl hist msg)
  (begin
;    (if (string=? msg "")
       (shdot)
;      (format #t "~a " msg)
;    )
    (if (in-list? f hist)
        f
        (let ((rf (replace-simple-forms f t vl))
              (hist1 (cons f hist))
             )
          (if (not (in-list? rf hist1))
              (replace-form rf t vl hist1 "")
              rf
          )
        )
    )
  )
))

(define comp-numbers1 (lambda (f)
  (letrec* ((c-n (lambda (fp t)
                   (let ((vrs (form-vars fp t (list) (make-hash-table))))
                        ;(if (and (list? fp) (not (null? fp))
                        ;         ;;(eq? (car fp) 'sign)
                        ;         )
                        ;    (format #t "SIGN: ~a\nVRS: ~a\n" (car fp) vrs)
                        ;)
                        (if (unspecified? vrs)
                            fp
                            (cond ((eq? (car fp) '**)
                                   (apply expt (cdr fp))
                                  )
                                  ((eq? (car fp) 'neg)
                                   (* -1.0 (cadr fp))
                                  )
                                  ((eq? (car fp) 'inv)
                                   (/ 1.0 (cadr fp))
                                  )
                                  ((eq? (car fp) 'sign)
                                   (cond ((> (cadr fp) 0.0) 1.0)
                                         ((< (cadr fp) 0.0) -1.0)
                                         (#t 0.0)
                                   )
                                  )
                                  (#t (primitive-eval fp))
                            )
                        )
                   )
            ))
            (c-n-l (lambda (fp l)
              (if (null? l)
                  fp
                  (c-n-l (c-n fp (car l)) (cdr l))
              )
            ))
            (c-n-f (lambda (fp)
              (if (and (list? fp))
                  (c-n-l (map c-n-f fp) math-num-forms)
                  fp
              )
            ))
          )
    (c-n-f f)
  )
))

(define comp-numbers (lambda (v hist)
  (shdot)
  (let ((v1 (comp-numbers1 v))
        (hist1 (cons v hist))
       )
    (if (not (in-list? v1 hist))
        (comp-numbers v1 hist1)
        v1
    )
  )
))

(define expand-form (lambda (v)
    (comp-numbers
      (replace-form (expand-to-wires v) 
                     expand-templates
                    (list) (list) 
                    "Expanding:"
      )
      (list)
    )
))

(define sort-form2 (lambda (v vl)
  (replace-form v 
                sort-templates
                vl
                (list)
                "Sorting:"
  )
))

(define sort-form (lambda (v vl)
  (if (< (length vl) 2)
      (sort-form2 v vl)
      (sort-form2 (sort-form v (cdr vl)) vl)
  )
))

(define easy-mul (lambda (v)
  (let* ((vl (filter (lambda (x) (not (and (number? x) (= x 1)))) (cdr v)))
         (zl (filter (lambda (x) (and (number? x) (= x 0))) vl))
         (nnl (filter (lambda (x) (not (number? x))) vl))
         (nl  (filter (lambda (x) (number? x)) vl))
         (nm  (apply * nl))
        )
    (if (not (null? zl))
        0.0
        (if (null? vl)
            1.0
            (if (= (length vl) 1)
                (car vl)
                (if (null? nnl)
                    nm
                    (cons (car v) (cons nm nnl))
                )
            )
        )
    )
  )
))

(define easy-plus (lambda (v)
  (let* ((vl (filter (lambda (x) (not (and (number? x) (= x 0)))) (cdr v)))
         (nnl (filter (lambda (x) (not (number? x))) vl))
         (nl  (filter (lambda (x) (number? x)) vl))
         (np  (apply + nl))
        )
    (if (null? vl)
        0.0
        (if (= (length vl) 1)
            (car vl)
            (if (null? nnl)
                np
                (cons (car v) (cons np nnl))
            )
        )
    )
  )
))

(define easy-minus (lambda (v)
  (let* ((vl  (easy-plus (cons '+ (cddr v))))
         (vl1 (if (and (list? vl) (not (null? vl)) (eq? (car vl) '+))
                  (cdr vl)
                  (list vl)
              )
         )
         (nnl (filter (lambda (x) (not (number? x))) vl1))
         (nl  (filter (lambda (x) (number? x)) vl1))
         (np  (apply + nl))
        )
    (if (and (number? (cadr v)) (null? nnl))
      (apply - (cons (cadr v) nl))
      (cons (car v) (cons (cadr v) (cons np nnl)))
    )
  )
))

(define easy-div (lambda (v)
  (if (and (number? (cadr v)) (= (cadr v) 0))
      0.0
      (if (and (number? (caddr v)) (= (caddr v) 1) (= (length v) 3))
          (cadr v)
          v
      )
  )
))

(define easy-comp (lambda (v)
  (if (and (list? v) (not (null? v)))
      (let* ((cv (car v))
             (v1 (cons cv (map easy-comp (cdr v))))
            )
        (cond ((eq? cv '+) (easy-plus v1))
              ((eq? cv '-) (easy-minus v1))
              ((eq? cv '*) (easy-mul v1))
              ((eq? cv '/) (easy-div v1))
              (#t v1)
        )
      )
      v
  )
))

(define easy-simplify-form2 (lambda (v hist)
  (shdot)
  (let ((ef (replace-form
              (comp-numbers v (list))
                simp1-templates
                (list)
                (list) 
                "Simp1:"
              )
        )
        (hist1 (cons v hist))
       )
    (if (not (in-list? ef hist1))
        (easy-simplify-form2 ef hist1)
        ef
    )
  )
))

(define easy-simplify-form (lambda (v)
  (norm-form
    (easy-simplify-form2 (expand-to-wires v) (list))
  )
))



(define get-all-vars1 (lambda (v)
  (cond ((symbol? v) (list v))
        ((procedure? v) (list v))
        ((and (list? v) (> (length v) 1))
         (append (get-all-vars1 (cadr v))
                 (get-all-vars1 (cons 1 (cddr v)))
         )
        )
        (#t (list))
  )
))

(define count-var (lambda (v l)
  (cond ((and (or (symbol? l) (procedure? l)) (equal? v l)) 1)
        ((and (list? l) (not (null? l)))
         (+ (count-var v (car l)) (count-var v (cdr l)))
        )
        (#t 0)
  )
))

(define get-all-vars (lambda (v)
  (let* ((vrs (remove-duplicates (get-all-vars1 v)))
         (cvrs (zip (map (lambda (x) (count-var x v)) vrs) vrs))
         (svrs (sort cvrs (lambda (x y) (< (car x) (car y)))))
        )
    (map cadr svrs)
  )
))

(define upd-nums (lambda (v)
  (cond ((procedure? v) v)
        ((number? v) (string->number (number->string v)))
        ((list? v) (map upd-nums v))
        (#t v)
  )
))

(define simplify2-form (lambda (v vrs hist)
  (let* ((rf1 (replace-form (comp-numbers (sort-form v vrs) (list))
                            (append
                              simplify1-templates
                              simplify2-templates
                            )
                            (list)
                            (list)
                            "Simplify2:"
               )
         )
         (rf rf1);(upd-nums rf1))
         (hist1 (cons v hist))
        )
    (if (not (in-list? rf hist1))
        (simplify2-form rf vrs hist1)
        rf
    )
  )
))

(define simplify3-form (lambda (v vrs)
  (if (null? vrs)
      v
      (simplify3-form (simplify2-form v vrs (list))
                      (reverse (cdr (reverse vrs)))
      )
  )
))

(define l->s (lambda (l)
  (cond ((string? l) (format #f "\"~a\"" l))
        ((list? l)
         (format #f "(~a)" (string-join (map l->s l) " "))
        )
        (#t (format #f "~a" l))
  )
))

(define simplify4-form (lambda (v vrs hist)
;  (format #t "\nSF4: ~a\n" (l->s v))
;  (format #t "\nVRS: ~a\n" (l->s vrs))
  (let* ((rf (simplify3-form v vrs))
         (erf rf);(expand-form rf))
         (hist1 (cons v hist))
        )
;    (format #t "RES: ~a\n" (l->s rf))
    (if (in-list? erf hist1)
        rf
        (simplify4-form erf vrs hist1)
    )
  )
))

(define simplify-form (lambda (v)
  (letrec* ((v1 (simp-ex v))
            (res (simplify4-form (expand-form v1) (get-all-vars v1) (list)))
            (p-vrs (lambda (v)
              (cond ((is-var? v) (var-name v))
                    ((list? v)
                     (map p-vrs v)
                    )
                    (#t v)
              )
            ))
           )
;    (format #t "\n\nSF: ~a\n" (p-vrs v))
;    (format #t "RES: ~a\n\n" (p-vrs res))
    (norm-form res)
  )
))

(def-cache norm-form (lambda (v)
  (replace-form 
    (replace-form v norm1-templates (list) (list) "Norm1:")
    norm2-templates (list) (list) "Norm2:"
  )
))

;; Solve e(v)=0 for v

(define solve-eq (lambda (e v)
  (letrec* ((e1 (sort-form (expand-form e) (list v)))
            ;;(xxy (format #t "\n\nE: ~a\n\n" e))
            ;;(xxx (format #t "\n\nE1[~a]: ~a\n\n" v e1))
            (try-form (lambda (t)
              (let ((vrs (form-vars e1 (car t) (list v) (make-hash-table))))
                ;;(format #t "SVRS: ~a\n" vrs)
                (if (not (unspecified? vrs))
                    (let ((res (subst-vars (cadr t) vrs)))
                         (if (not (cont-one res (list v)))
                             res
                         )
                    )
                )
              )
            ))
            (try-forms (lambda (l)
              (if (not (null? l))
                  (let ((res (try-form (car l))))
                       (if (not (unspecified? res))
                           res
                           (try-forms (cdr l))
                       )
                  )
              )
            ))
            (res (try-forms solve-templates))
           )
    ;;(if (unspecified? res)
    ;;    res
        (simp-ex res)
    ;)
  )
))

(define has-func (lambda (v f)
  (letrec ((h-f-l (lambda (l)
             (and (not (null? l))
                  (or (has-func (car l) f)
                      (h-f-l (cdr l))
                  )
             )
          )))
    (and (list? v)
         (not (null? v))
         (or (eq? (car v) f)
             (h-f-l (cdr v))
         )
    )
  )
))

;; df(x)/dx

(def-cache diff-eq (lambda (f x)
  (let* ((res (replace-form (list 'diff (expand-form f))
                            diff-templates
                            (list x)
                            (list)
                            "Diff:"
              )
         )
         ;;(sres (simp-ex res))
        )
    (if (has-func res 'diff)
        (int-err (format #f "Cannot diff (~a): ~a" x res))
        res
    )
  )
))

(define extract-consts (lambda (e vrs tvs)
  (begin
  (cond ((number? e) (list (list) e))
        ((not (cont-one e vrs))
         (list (list (list (string->symbol tvs) 
                           (simp-ex 
                             e
                           )
                     )
               )
               (string->symbol tvs)
         )
        )
        ((and (list? e) (= (length e) 2))
         (let ((ec1 (extract-consts (cadr e)  vrs (string-append tvs "0")))
              )
              (list (car ec1)
                    (list (car e) (cadr ec1))
              )
         )
        )
        ((and (list? e) (= (length e) 3))
         (let ((ec1 (extract-consts (cadr e)  vrs (string-append tvs "0")))
               (ec2 (extract-consts (caddr e) vrs (string-append tvs "1")))
              )
              (list (append (car ec1) (car ec2))
                    (list (car e) (cadr ec1) (cadr ec2))
              )
         )
        )
        (#t (list (list) e))
        
  )
  )
))

;;; Solve system Y[X]=0 for X

(define solve-sys (lambda (y x)
  (letrec* ((try-all (lambda (yl xl x0)
                       (if (and (not (null? xl))
                                (not (null? yl))
                                (cont-one (car yl) (list (car xl)))
                           )
                         (let ((sl (solve-eq ;;(simp-ex
                                              (car yl)
                                             ;;)
                                             (car xl)
                                   )
                               )
                              )
                              (if (unspecified? sl)
                                  (try-all yl (cdr xl) x0)
                                  (list (car xl) sl (car yl))
                              )
                         )
                         (if (not (null? yl))
                             (try-all (cdr yl) x0 x0)
                         )
                       )
                     )
            )
            (solve-all (lambda (yl xl rl)
                         (if (null? xl)
                             rl
                             (let* ((sl (try-all yl xl xl))
                                    (ssl (if (not (unspecified? sl))
                                             (list
                                               (car sl)
                                               (simp-ex 
                                                 (cadr sl)
                                               )
                                               (caddr sl)
                                             )
                                         )
                                    )
                                   )
                               ;;(format #t "SL: ~a\n" sl)
                               ;;(format #t "SSL: ~a\n" ssl)
                               (if (unspecified? ssl)
                                   (int-err "Cannot solve system!")
                                   (let ((ht (make-hash-table)))
                                     (hash-set! ht (car ssl) (cadr ssl))
                                     (hash-for-each
                                       (lambda (f v)
                                         (hash-set! ht f (subst-vars v ht))
                                       )
                                       rl
                                     )
                                     (solve-all
                                       (map (lambda (f) 
                                              (subst-vars f ht)
                                            )
                                            (delete (caddr ssl) yl)
                                       )
                                       (delete (car ssl) xl)
                                       ht
                                     )
                                   )
                               )
                             )
                         )
                       )
            )
            (e-c-l (lambda (l)
                     (if (not (null? l))
                         (let ((cl1 (extract-consts 
                                      (car l)
                                      (arr-elements x)
                                      (format #f "~~CNS~a_" (length l))
                                    )
                               )
                               (cl2 (e-c-l (cdr l)))
                              )
                           (list (append (car cl1) (car cl2))
                                 (cons (cadr cl1) (cadr cl2))
                           )
                         )
                         '(() ())
                     )
                   )
            )
            (res (solve-all 
                   (arr-elements y)
                   (arr-elements x)
                   (make-hash-table)
                 )
            )
            (h-res (make-hash-table))
           )
    (hash-for-each
      (lambda (f v)
        (hash-set! h-res f (simp-ex (norm-form (subst-vars v res))))
      )
      res
    )
    h-res
  )
))

(define temps (lambda (l) (map lamcon (iota l))))


(define jakobi-matrix (lambda (f q vrs)
  (letrec* ((tmp-vars (temps (cadr (math-array-dims f))))
            (qres (solve-sys (apply math-cvector 
                                    (map 
                                      (lambda (x) `(- ,(cadr x) ,(car x)))
                                      (zip tmp-vars (arr-elements q))
                                    )
                             )
                             (apply math-cvector vrs)
                  )
            )
            (fq (subst-vars (arr-elements f) qres))
            (rdims (list
                     (cadr (math-array-dims f))
                     (cadr (math-array-dims q))
                   )
            )
            (max-i (* (car rdims) (cadr rdims)))
            (make-el (lambda (i)
              (let ((ind (arr-n->ind i rdims)))
                (if (< i max-i)
                    (let* ((df (diff-eq 
                                 (list-ref fq (cadr ind))
                                 (list-ref tmp-vars (car ind))
                               )
                           )
                           (thl (make-hash-table))
                           (rs (begin
                                 (for-each
                                   (lambda (x)
                                     (hash-set! thl (car x) (cadr x))
                                   )
                                   (zip tmp-vars (arr-elements q))
                                 )
                                 (subst-vars df thl)
                               )
                           )
                           (rss (simp-ex
                                  rs
                                )
                           )
                          )
                      ;(format #t "J[~a:~a]= ~a\n" 
                      ;  (car ind) (cadr ind) 
                      ;  rss
                      ;  ;;(infix-out rss)
                      ;)
                      (cons rss (make-el (+ i 1)))
                    )
                    (list)
                )
              )
            ))
           )
;    (format #t "F1: ~a\n" (infix-out (car fq)))
;    (format #t "F2: ~a\n" (infix-out (cadr fq)))
;    (format #t "F3: ~a\n" (infix-out (caddr fq)))
    (math-array rdims (make-el 0))
  )
))


(define get-right-ev (lambda (a l tmpl)
  (let* ((dims (math-array-dims a))

         (tmp-r (temps (apply * dims)))
         (rt (math-array dims tmp-r))
         (rrt (apply math-cvector
                      (map
                        car
                        (filter (lambda (x) (eq? (cadr x) '?))
                                (zip (arr-elements rt) tmpl)
                        )
                      )
              )
         )
         (srt (filter (lambda (x) (not (eq? (cadr x) '?)))
                      (zip (arr-elements rt) tmpl)
              )
         )
         
         (hsrt (make-hash-table))
         
         (lam (diag l))
         (ar (mul2 a rt))
         (rl (mul2 rt lam))
         (eqs (math-array
                dims 
                (map (lambda (x) (subst-vars x srt))
                     (arr-elements (add ar (mul2 -1.0 rl)))
                )
              )
         )
         (sl (solve-sys eqs rrt))
        )
    (for-each
      (lambda (x)
        (hash-set! hsrt (car x) (cadr x))
      )
      srt
    )
    (hash-for-each
      (lambda (x v)
        (hash-set! hsrt x (simp-ex (subst-vars v hsrt)))
      )
      sl
    )

    (if (not (unspecified? sl))
      (math-array dims (subst-vars (arr-elements rt) hsrt))
      (int-err "Cannot find right eigenvectors")
    )
  )
))

;(define get-left-ev (lambda (a l tmpl)
;  (let* ((dims (math-array-dims a))
;
;         (tmp-l (temps (apply * dims)))
;         (lt (math-array dims tmp-l))
;         (llt (apply math-cvector
;                      (map
;                        car
;                        (filter (lambda (x) (eq? (cadr x) '?))
;                                (zip (arr-elements lt) tmpl)
;                        )
;                      )
;              )
;         )
;         (slt (filter (lambda (x) (not (eq? (cadr x) '?)))
;                      (zip (arr-elements lt) tmpl)
;              )
;         )
;         
;         (lam (diag l))
;         (la (mul2 lt a))
;         (ll (mul2 lam lt))
;         (eqs (math-array
;                dims 
;                (map (lambda (x) (subst-vars x slt))
;                     (arr-elements (add la (mul2 -1.0 ll)))
;                )
;              )
;         )
;         (sl (solve-sys eqs (math-array dims (reverse (arr-elements llt)))))
;        )
;    (if (not (unspecified? sl))
;      (math-array dims 
;                  (subst-vars
;                    (arr-elements lt)
;                    (append
;                      (map 
;                        (lambda (x)
;                          (list (car x)
;                                (simp-ex 
;                                  (subst-vars (cadr x) slt)
;                                )
;                          )
;                        )
;                        sl
;                      )
;                      slt
;                    )
;                  )
;      )
;      (int-err "Cannot find left eigenvectors")
;    )
;  )
;))

;(define inverse-matrix (lambda (x)
;  (let* ((dims (math-array-dims x))
;  
;         (ones (apply math-cvector (map (lambda (x) 1.0) (iota (car dims)))))
;         (I (diag ones))
;
;         ;(tmp-x (temps (apply * dims)))
;         ;(tx (math-array dims tmp-x))
;         ;(ts (zip tmp-x (arr-elements x)))
;         
;         (tmp-y (temps (apply * dims)))
;         (ty (math-array dims tmp-y))
;
;         (eqs (math-array dims ;;(map simp-ex 
;                      (arr-elements (add (mul2 ty x) (mul2 -1.0 I)))
;                               ;;)
;             ))
;         (sl (solve-sys eqs ty))
;        )
;    (if (not (unspecified? sl))
;      (math-array dims 
;                  (map simp-ex
;                      ;(subst-vars
;                        (subst-vars
;                          (arr-elements ty)
;                          sl
;                        )
;                      ;  ts
;                      ;)
;                  )
;      )
;      (int-err "Cannot find inverse matrix")
;    )
;  )
;  
;))

(define sym-det (lambda (ar)
  (if (apply = (math-array-dims ar))
    (simp-ex
      (letrec*
        ((els (arr-elements ar))
         (dims (math-array-dims ar))
         (rank (car dims))
         (inds (all-inds (list 0) (map (lambda (x) rank) (iota1 rank))))
         (lp-v (lambda (i)
           (if (null? i)
               0.0
               (let ((lc (levi-civita (car i))))
                 (if (not (= lc 0))
                    `(+ ,(cons '*
                           (cons lc
                                 (map
                                   (lambda (x)
                                     (list-ref
                                       els
                                       (arr-ind->n x dims)
                                     )
                                   )
                                   (zip (car i) (iota1 rank))
                                 )
                           )
                         )
                        ,(lp-v (cdr i))
                     )
                     (lp-v (cdr i))
                 )
               )
           )
         ))
        )
        (if (= rank 1)
            (car els)
            (lp-v inds)
        )
      )
    )
    (int-err "Det: matrix must be square!")
  )
))

(define sym-vec-mul (lambda (v1a v2a)
  (letrec*
    ((n-v (lambda (v)
       (let* ((d (math-array-dims v))
              (md (apply max d))
              (mdi (list-index (lambda (x) (> x 1)) d))
              (mdi1 (if mdi mdi (- (length d) 1)))
             )
         (if (> (length (filter (lambda (x) (> x 1)) d)) 1)
             (int-err "Vec-mul: not a vector")
             (if (< md 3)
                 (math-array
                   (replace-el d mdi1 3)
                   (append
                     (arr-elements v)
                     (map 
                       (lambda (x) 0)
                       (iota1 (- 3 md))
                     )
                   )
                 )
                 v
             )
         )
       )
     ))
     (v1 (n-v v1a))
     (v2 (n-v v2a))
     (dims  (math-array-dims v1))
     (dims2 (math-array-dims v2))
     (els1 (arr-elements v1))
     (els2 (arr-elements v2))
     (rank (apply max dims))
     (inds (all-inds (list 0) (list rank rank)))
     (l-p (lambda (i)
       (if (< i rank)
           (cons
             (letrec
               ((lp-ind (lambda (ind)
                  (if (null? ind)
                      0.0
                     `(+ (* ,(levi-civita (cons i (car ind)))
                            ,(list-ref els1 (car (car ind)))
                            ,(list-ref els2 (cadr (car ind)))
                         )
                        ,(lp-ind (cdr ind))
                      )
                  )
                ))
               )
               (lp-ind inds)
             )
             (l-p (+ i 1))
           )
           (list)
       )
     ))
    )
    (if (equal? dims dims2)
        (math-array dims (map simp-ex (l-p 0)))
        (int-err "Vec-mul: invalid parameters")
    )
  )
))

(define sym-vec-mod (lambda (v)
  (simp-ex
    (list 'sqrt (cons '+ (map (lambda (x) `(** ,x 2.0)) (arr-elements v))))
  )
))

;(define grad (lambda (f bas vrs)  ;; untested
;  (arr-row (jakobi-matrix (math-cvector f 0 0) bas vrs) 0)
;))

(define get-nums (lambda (l)
  (filter number? l)
))

(define get-non-nums (lambda (l)
  (filter (lambda (x) (not (number? x))) l)
))

(define join-same-funcs (lambda (f)
  (lambda (v)
    (if (or (number? v) (procedure? v) (is-ptr? v) (symbol? v))
        (list v)
        (if (and (list? v) (> (length v) 1) (eq? (car v) f))
            (apply append (map (join-same-funcs f) (cdr v)))
            (list v)
        )
    )
  )
))

(def-cache simp-ex2 (lambda (ex)
  (cond ((is-ptr? ex) ex)
        ((is-var? ex)
         (if (var-has-ind? ex)
             (let* ((vi (var-ind ex))
                    (vi2 (map simp-ex2 vi))
                   )
               (set-ind (del-ind ex) vi2)
             )
             ex
         )
        )
        ((or (is-op? ex) (is-fn? ex))
         (let ((exv (map simp-ex2 (cdr ex))))
           (cond ((eq? (car ex) '+)
                  (let* ((exv1 (apply append (map (join-same-funcs '+) exv)))
                         (nm (apply + (get-nums exv1)))
                         (nnm (get-non-nums exv1))
                        )
                    (if (null? nnm)
                        nm
                        (if (= nm 0.0)
                            (if (= (length nnm) 1)
                                (car nnm)
                                (cons (car ex) nnm)
                            )
                            (cons (car ex) (append nnm (list nm)))
                        )
                    )
                  )
                 )
                 ((eq? (car ex) '-)
                  (if (= (length exv) 2)
                    (if (and (number? (car exv)) (number? (cadr exv)))
                        (apply - exv)
                        (if (and (number? (cadr exv)) (= (cadr exv) 0))
                            (car exv)
                            (cons (car ex) exv)
                        )
                    )
                    (let* ((exv1 (cons
                                   (car exv)
                                   (apply append (map (join-same-funcs '+) 
                                                     (cdr exv)))
                                 )
                           )
                           ;;(exv1 exv)
                           (nm (apply + (get-nums (cdr exv1))))
                           (nnm (get-non-nums (cdr exv1)))
                           (rs1 (if (= nm 0.0)
                                    (cons (car exv) nnm)
                                    (cons (car exv) (append nnm (list nm)))
                                )
                           )
                           (nnm2 (get-non-nums rs1))
                          )
                      (if (null? nnm2)
                          (if (= (length rs1) 1)
                              (car rs1)
                              (apply - rs1)
                          )
                          (if (= (length rs1) 1)
                              (car rs1)
                              (cons (car ex) rs1)
                          )
                      )
                    )
                  )
                 )
                 ((eq? (car ex) '*)
                  (let* ((exv1 (apply append (map (join-same-funcs '*) exv)))
                         (nm (apply * (get-nums exv1)))
                         (nnm (get-non-nums exv1))
                        )
                    (if (null? nnm)
                        nm
                        (if (= nm 0.0)
                            0.0
                            (if (= nm 1.0)
                                (if (= (length nnm) 1)
                                    (car nnm)
                                    (cons (car ex) nnm)
                                )
                                (cons (car ex) (append nnm (list nm)))
                            )
                        )
                    )
                  )
                 )
                 ((eq? (car ex) '/)
                  (let ((nms (get-nums exv))
                        (nnm (get-non-nums exv))
                       )
                    (if (null? nnm)
                        (if (and-map exact? nms)
                            (cons (car ex) exv)
                            (apply / nms)
                        )
                        (if (and (number? (car exv)) (= (car exv) 0))
                            0.0
                            (if (and (= (length exv) 2)
                                     (number? (cadr exv))
                                     (= (cadr exv) 1)
                                )
                                (car exv)
                                (cons (car ex) exv)
                            )
                        )
                    )
                  )
                 )
                 ((eq? (car ex) '**)
                  (if (and (= (length exv) 2) (number? (cadr exv)))
                      (if (= (cadr exv) 0)
                          1.0
                          (if (= (cadr exv) 1)
                              (car exv)
                              (if (number? (car exv))
                                  (expt (car exv) (cadr exv))
                                  (cons (car ex) exv)
                              )
                          )
                      )
                      (cons (car ex) exv)
                  )
                 )
                 ((eq? (car ex) 'abs)
                  (if (and (= (length exv) 1) (number? (car exv)))
                      (abs (car exv))
                      (cons (car ex) exv)
                  )
                 )
                 ((eq? (car ex) 'sign)
                  (if (and (= (length exv) 1) (number? (car exv)))
                      (if (< (car exv) 0) -1.0 1.0)
                      (cons (car ex) exv)
                  )
                 )
                 ((eq? (car ex) 'max)
                  (let* ((exv1 (apply append (map (join-same-funcs 'max) exv)))
                         (nms (get-nums exv1))
                         (nm (if (null? nms) #f (apply max nms)))
                         (nnm (remove-duplicates (get-non-nums exv1)))
                        )
                    (if (null? nnm)
                        nm
                        (cons (car ex) (append nnm (if nm (list nm) '())))
                    )
                  )
                 )
                 ((eq? (car ex) 'min)
                  (let* ((exv1 (apply append (map (join-same-funcs 'min) exv)))
                         (nms (get-nums exv1))
                         (nm (if (null? nms) #f (apply min nms)))
                         (nnm (remove-duplicates (get-non-nums exv1)))
                        )
                    (if (null? nnm)
                        nm
                        (cons (car ex) (append nnm (if nm (list nm) '())))
                    )
                  )
                 )
                 ((eq? (car ex) 'sqrt)
                  (if (and (= (length exv) 1) (number? (car exv)))
                      (sqrt (car exv))
                      (cons (car ex) exv)
                  )
                 )
                 (#t (cons (car ex) exv))
           )
         )
        )
        (#t ex)
  )
))

(def-cache mk-mul-sum (lambda (ex)
  (cond ((or (is-fn? ex) (is-op? ex))
         (let ((exv (map mk-mul-sum (cdr ex))))
           (cond ((eq? (car ex) '-)
                  (cons '+
                        (cons (car exv)
                              (map (lambda (x) `(* ,x -1))
                                   (cdr exv)
                              )
                        )
                  )
                 )
                 ((eq? (car ex) '/)
                  (cons '*
                        (cons (car exv)
                              (map (lambda (x) `(/ 1 ,x))
                                   (cdr exv)
                              )
                        )
                  )
                 )
                 (#t (cons (car ex) exv))
           )
         )
        )
        ((and (is-var? ex) (var-has-ind? ex))
         (set-ind (del-ind ex) (map mk-mul-sum (var-ind ex)))
        )
        (#t ex)
  )
))

(def-cache mk-minus-div (lambda (ex)
  (cond ((or (is-op? ex) (is-fn? ex))
         (let ((exv (map mk-minus-div (cdr ex))))
           (cond ((and (eq? (car ex) '+)
                       (= (length exv) 2)
                       (number? (cadr exv))
                       (< (cadr exv) 0)
                  )
                 `(- ,(car exv) ,(* (cadr exv) -1))
                 )
                 ((and (eq? (car ex) '+)
                       (> (length exv) 1)
                  )
                  (letrec* ((neg-pred (lambda (x)
                                        (and (or (is-op? x) (is-fn? x))
                                             (> (length x) 2)
                                             (eq? (car x) '*)
                                             (number? (car (reverse x)))
                                             (= (car (reverse x)) -1)
                                        )
                                      )
                            )
                            (nneg-pred (lambda (x) (not (neg-pred x))))
                            (negs (filter neg-pred exv))
                            (nnegs (filter nneg-pred exv))
                            (chp (lambda (x)
                               (if (> (length x) 3)
                                   (reverse (cdr (reverse x)))
                                   (cadr x)
                               )
                            ))
                           )
                    (if (null? negs)
                        (cons '+ exv)
                        (if (null? nnegs)
                           `(* ,(cons '+ (map chp negs))
                               ,(car (reverse (car negs)))
                            )
                           `(- ,(if (= (length nnegs) 1)
                                    (car nnegs)
                                    (cons '+ nnegs)
                                )
                               ,(if (= (length negs) 1)
                                    (chp (car negs))
                                    (cons '+ (map chp negs))
                                )
                            )
                        )
                    )
                  )
                 )
                 ((and (eq? (car ex) '*)
                       (> (length exv) 1)
                  )
                  (letrec* ((inv-pred (lambda (x)
                                        (and (or (is-op? x) (is-fn? x))
                                             (> (length x) 2)
                                             (eq? (car x) '/)
                                             (number? (cadr x))
                                             (= (cadr x) 1)
                                        )
                                      )
                            )
                            (ninv-pred (lambda (x) (not (inv-pred x))))
                            (invs (filter inv-pred exv))
                            (ninvs (filter ninv-pred exv))
                           )
                    (if (null? invs)
                        (cons '* exv)
                        (if (null? ninvs)
                           `(/ ,(cadr (car invs))
                               ,(cons '* (map caddr invs))
                            )
                           `(/ ,(if (= (length ninvs) 1)
                                    (car ninvs)
                                    (cons '* ninvs)
                                )
                               ,(if (= (length invs) 1)
                                    (caddr (car invs))
                                    (cons '* (map caddr invs))
                                )
                            )
                        )
                    )
                  )
                 )
                 (#t (cons (car ex) exv))
           )
         )
        )
        ((and (is-var? ex) (var-has-ind? ex))
         (set-ind (del-ind ex) (map mk-minus-div (var-ind ex)))
        )
        (#t ex)
  )
))

(def-cache simp-ex (lambda (ex)
  (mk-minus-div (simp-ex2 (mk-mul-sum ex)))
))

)
